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We study the relaxation of a non-equilibrium carrier dis- 
tribution under the influence of the electron-electron in- 
teraction in the presence of disorder. 
Based on the Anderson model, our Hamiltonian is com- 
posed from a single particle part including the disorder 
and a two-particle part accounting for the Coulomb in- 
teraction. We apply the equation-of-motion approach for 
the density matrix, which provides a fully microscopic 
description of the relaxation. 



Our results show that the nonequlibrium distribution in 
this closed and internally interacting system relaxes ex- 
ponentially fast during the initial dynamics. This fast re- 
laxation can be described by a phenomenological damp- 
ing rate. The total single particle energy decreases in 
the redistribution process, keeping the total energy of 
the system fixed. It turns out that the relaxation rate de- 
creases with increasing disorder. 
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1 Introduction Understanding the interplay of disor- 
der and interaction constitutes a challenging problem in 
condensed matter physics [ 1 ]. Even weak interactions cause 
a considerable amount of complexity [2]. One of the most 
intriguing phenomena, the phase relaxation due to inter- 
action, has already been treated in the presence of disor- 
der [3]. In many cases, however, either of the two major 
ingredients, the disorder and the inter-particle interaction 
are taken into account only perturbatively. In numerical 
simulations this type of approximation is not necessary. 
However, finite-size limitations do affect the results, even 
though further controlled simplifications may nevertheless 
be necessary to impose. 

In this paper we present a case study where we have 
investigated how the interplay of disorder and long-range 
interaction affects the relaxation of a non-equilibrium en- 
ergy distribution of electrons in one dimension. 

In this paper our aims are twofold. On the one hand, we 
want to determine whether the relaxation within a closed 
interacting system can be described by a phenomenologi- 
cal damping rate. On the other hand, we perform the first 
steps in the direction of constructing a theory of relaxation 
treating both disorder and interaction on an equal micro- 



scopic footing. Our approach is based on the equation of 
motion technique applied on elements of the density ma- 
trix. This technique has achieved a great success in nu- 
merous applications in semiconductor optics for systems 
with strong many -body correlations [4,5]. Our experience 
shows that the interplay of disorder and interaction does 
yield novel and interesting phenomena [6,7]. Here we pro- 
vide one more evidence that it can also be used to tackle 
challenging problems in the field of interacting disordered 
systems. 

The phase relaxation due to inter-particle interactions 
in disordered systems has already been investigated analyt- 
ically, for a recent review, see [3]. As a summary, the pres- 
ence of weak disorder accelerates the relaxation process, 
since the scattering probabilities increase due to the lack 
of momentum selection rules an ordered system would im- 
pose. The question of strong disorder, however, is not yet 
fully understood. A first step using numerical simulations 
was performed in [8], where a Boltzmann-type equation 
was solved with an initial condition of a non-equilibrium 
carrier distribution, which relaxes over the course of time 
towards a Fermi-Dirac-type distribution. While the analyt- 
ical prediction for weak disorder was recovered in the case 
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when short-ranged interactions were allowed between the 
particles, for both long ranged and short ranged interaction 
the strong disorder limit yielded a decrease in the relax- 
ation rate. The latter behavior was attributed to the small 
size of the localization volume of the one particle-states, 
resulting in smaller scattering probabilities. 

Since the previous work of I. Varga et al. [8] suffered 
from several simplifying assumptions, in the present work 
we perform computations along similar ideas, but using a 
much better founded formalism, namely the equation of 
motion of the density matrix [5, 6, 7], 

2 Model 

2.1 The Hamiltonian Our model is based on a one 
dimensional tight-binding model with energetic disorder, 
represented by the energies of the sites denoted by £j for 
the j*^ site [5,7]. The model also includes nearest neighbor 
hopping matrix elements Jij ~ J between sites i and j. 
Using a site basis where \i) denotes the basis function at 
site i it is straightforward to formulate the non-interacting 
part of the Hamiltonian: 



He 



<ij> 



where < ij > corresponds to a sum limited to the nearest 
neighbors. 

The term responsible for the two-body interaction reads 

[5] 



He 



\j) 01 (*l 



(2) 



where only monopole-monopole terms have been included. 
In Eq. (2), the sum represents the repulsive Coulomb inter- 
action among charge carriers within a single band. The reg- 
ularized Coulomb matrix element considered here is given 
by 

V.,= .. ^° , (3) 
\i- j\a + ao 

with a positive constant Uq. Here a = 5 nm is the site sep- 
aration [5] and the term oq = 0.5 a removes the unphysical 
singularity of the lowest excitonic bound state arising from 
the restriction to monopole-monopole terms in one dimen- 
sion [5, 7]. 

2.2 Single Particle Basis In order to derive the equa- 
tions of motion it is useful to write the Hamiltonian in 
second-quantized formalism. As we are interested in the 
dynamics of single particles, we transform the Hamiltonian 
into the eigenbasis of single particles via diagonalizing Hq 
in Eq. (1) 

^o'^^E^"^"^"- (4) 

a 

After the transformations the Hamiltonian containing the 
Coulomb interaction reads as 



E 



V 



ct/3 C.^,Cp,CpCa 



(5) 



with the transformed matrix elements 



Here ^^(i) is the projection of the electronic eigenstate 
with quantum number a onto site j. 

3 The equation of motion The electron population 
in an eigenstate a is given by the expectation value of the 
number operator: 



The off-diagonal extensions 



(7) 



(8) 



are called coherences [5] and play an important role in the 
present calculation. The dynamics of all these matrix ele- 
ments of the density matrix are coupled via the Coulomb 
interaction. Thus one has to calculate the time evolution of 
the full density matrix, which leads to the following equa- 
tion of motion: 



i'h^^ni2 = (£2 - ei) "-12 



17/32 



(9) 



with the correlated two-particle quantity being 

Ca/375 = (c)^c^pC^cs^ , (10) 
and the symmetrized Coulomb matrix element defined as 



a' 13' 



a/3 



2 V'-P 



v: 



a'0' 



pa 



(11) 



Note that in Eq. (9) the indices 1 , 2 denote general single- 
particle states. They were introduced only to help the reader 
to separate the summing and non-summing indices within 
the equation. Due to the coupling of n.ap to Cap-^s in prin- 
ciple the equation of motion for the latter quantity should 
be derived as well, and the two coupled equations should 
be solved together. It is, however, possible to separate the 
contribution of uncorrelated particles and pure two-particle 
quantum correlation with the help of a factorization scheme, 
which yields [4] 

Cap-iS = UaSnpj — Ua^niiS — ACapjS ■ (12) 

Assuming the contribution of the pure two-particle part to 
be negligible in comparison with the uncorrelated terms, 
i.e., 

AC^p^s = 0, (13) 
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we end up with the equation of motion 

iJi^rii2 = (62 - ei) ni2 

af3^ 



(14) 



for the populations and the coupled off-diagonal elements. 

4 Results In order to solve Eq. (14) using standard 
numerical integration techniques, we start with an initial 
population distribution 



1 



(15) 



which is centered around the energy Ec with a half width 
at full maximum w. Z is fixed by the total number of elec- 
trons "-a ~ ^/2, where N is the number of sites, cho- 
sen to be = 20. The width of the distribution is chosen 
to be ui = J in order to avoid unphysical cases of too large 
values of n^. 

In order to quantify how much the system has relaxed 
we use two physical quantities. One of them is the total 
single particle energy 



E 



(16) 



which evolves in time as ria changes. The other quantity is 
a measure which tells us how far the system is from equi- 
librium. It is defined by the root mean square of standard 
deviation of a nonlinear fit to a Fermi-Dirac distribution: 



1/2 



(17) 



In general Eq. (14) can be solved only numerically for a 
particular realization of the disorder. Thus an averaging 
over numerous realizations has to be performed in order 
to eliminate the dependence on the specific disordered po- 
tential landscape. Applying this to the total energy yields 
the results shown in Fig. 1. The total single particle en- 
ergy decreases from its initial value with a rate depending 
on the mean level spacing of single-particle states A = 
(W+B) /N/fi (h = 0.658 meVps) that results in a roughly 
universal behavior when time is rescaled accordingly. 

The configuration average is possible only for quan- 
tities like the single particle energy, E^^ or the standard 
deviation, a. On the other hand, nap corresponds to the 
eigenstates a and /3 and the related eigenenergies, which 
are strongly dependent on the particular disorder realiza- 
tion. I.e., even the same a and /? refer to different eigen- 
states, which leads to non-comparable for different re- 
alizations. 
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Figure 1 Total single particle energy as a function of time. 
The inset shows the raw data, in the main panel time is 
rescaled with the mean level spacing. 
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Figure 2 Error to a fit to a Fermi-Dirac distribution for a 
fixed interaction strength, Ua = B = 32 meV for several 
values of disorder strength. 



Now let us turn towards the evaluation of the relax- 
ation rate, F, which is defined via the assumption that (j{t) 
decreases exponentially in the short-time regime of the dy- 
namics, i.e., 

a{t) = aoe~" . (18) 



Indeed, our assumption is corroborated by the numerical 
simulations, shown in Fig. 2. The figure shows the fall- 
off of the deviation cr as a function of time on a semi-log 
plot for several values of disorder in units of the bare band- 
width, B = 4 J. [As we have fixed the value of J = 8 meV 
the bare bandwidth is S = 32 meV. Hence the timestep for 
the integration in Eq. (14) was chosen to he St I fs.] 

Thus it is reasonable to determine F for various disor- 
der and interaction strength values, both of which are mea- 
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Figure 3 Relaxation Rate as a function of disorder 
strength for various interaction strengths. Data after rescal- 
ing as {Uo/ BY, c w 1.6 is plotted in the main panel. The 
inset shows the raw data as a function of W/B. 



sured in units of the bare bandwidth, B — 4 J. The results 
of these calculations are shown in Fig. 3, where F is plot- 
ted as a function of disorder strength. The different lines 
belong to different interaction strengths. The rate F is also 
given in units of B. 

The inset in Fig. 3 contains the raw data which after a 
rescaling as [Uq/BY is shown in the main panel, where 
c « 1.6. A similar scaling of the relaxation rate was found 
in [8], however with an exponent closer to the value of 2, 
but more work needs to be done to understand the nature 
of these scaling exponents. Clearly in Fig. 2 the case with 
low interaction strength deviates from the major trend sub- 
stantially. 

5 Summary In the present work we have shown that 
the interplay of both the disorder and the long-range inter- 
action can be investigated on the same footing in order to 
study the relaxation process of an initially non-equilibrium 
one-particle occupation distribution towards an equilibrium 
one. The numerical solution of the equation of motion of 
the density matrix yields an exponential form for the devi- 
ation of the distribution from a Fermi-Dirac one. Addition- 
ally we found that the total single-particle energy gradually 
decreases as a function of time. This way the relaxation 
process in a completely closed system is possible via the 
rearrangement of the energy into other components of the 
otherwise constant total energy. A similar result has been 
found in [8], which, however, contained a phenomenolog- 
ical coupling to an environment. Morevoer, in the present 
calculation the off -diagonal matrix elements of the density 
matrix play an important role. 

The results indicate that further advances can be achieved 
by improving the present method. In a subsequent work the 



effect of higher order correlations will be investigated, as 
well as the case of short ranged interaction. 
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